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Abstract 

The parameterization of small-scale turbulent fluctuations in convective systems and in the pres- 
ence of strong stratification is a key issue for many applied problems in oceanography, atmospheric 
science and planetology. In the presence of stratification, one needs to cope with bulk turbulent 
fluctuations and with inversion regions, where temperature, density -or both- develop highly non- 
linear mean profiles due to the interactions between the turbulent boundary layer and the unmixed 
-stable- flow above/below it. We present a second order closure able to cope simultaneously with 
both bulk and boundary layer regions, and we test it against high-resolution state-of-the-art 2D 
numerical simulations in a convective and stratified belt for values of the Rayleigh number, up to 
Ra ~ 10 9 . Data are taken from a Rayleigh- Taylor system confined by the existence of an adiabatic 
gradient. 



Realistic parameterization of convective regions in presence of strong stratification is a 
problem of interest for the evolution of the convective boundary layer in fields as different 
as atmospheric science [T], stellar convection [2J [3] and oceanography |3H7]. The problem is 
intriguing also from a pure theoretical point of view: one would like to predict the interplay 
between buoyancy and turbulence at the edge between a convective region and a stable 
stratified volume above/below it. Due to turbulent activity, intermittent puffs of temper- 
ature, density and momentum tend to enter the stably stratified region, locally producing 
an inversion in the energy balance: kinetic energy is indeed lost, as potential energy is in- 
creased. Those turbulent puffs are the result of intense plumes traveling across the volume, 
often generated in the bottom layer, producing a "non-local" interaction between edge and 
bulk of the turbulent medium. The problem, already recognized in the early 50's [HJ E], 
is still a subject of intense research (see e.g. [5]). Recent studies have shown that, in the 
absence of strong stratification, mixing length theory based on Prandtl or Spiegel closure 
works very well in situations as different as for the case of Rayleigh- Taylor systems [TU] or in 
fluid mixing by gravity currents [TT]. Nevertheless, whenever strong stratification stops the 
evolution of the mixing profile, an overshooting region with temperature inversion develops 
and local closures fail. The problem arises from non-local effects caused by turbulent plumes 
or convective updraft. 

The situation can be visualized in the middle panel of Fig. [T] where a highly turbulent 
configuration is confined by two stably stratified media below and above it, at an effective 
Rayleigh number, Ra ~ 10 9 . The mean temperature profile is linear in the bulk and it 
develops two -symmetric- overshooting regions at the edge between the turbulent boundary 
layer and the fluid at rest with the heat flux inverting sign (see Fig. [2|. The configuration 
is obtained by evolving a Rayleigh- Taylor (RT) system in 2D with stratification [X2TTT7] . 
The evolution of the RT mixing layer, unbounded in the absence of stratification, is stopped 
because of the adiabatic gradient, i.e. when the mean profile of the potential temperature 
becomes flat. In this letter we propose a closure for the evolution of the mixing layer, able 
to capture both the initial transient (free of any stratification effect) and the late slowing 
down and stopping due to stratification effects. The main idea is to go beyond single point 
closure for the mean temperature evolution, keeping it exact and closing only for second order 
quantities: the Reynolds stresses, the heat flux and temperature fluctuations, see e.g. [3]. 
We test our closure against state-of-the-art numerical simulations at high resolution. We 
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FIG. 1: Right panel: initial vertical density and temperature profiles. Left and middle panels: 
Vorticity and temperature snapshots with weak and strong stratification, respectively. Notice 
the free rising plumes in the left panel typical of non-stratified RT system. In the middle panel 
stratification is stronger and turbulence is confined below and above the unmixed fluid at rest: at 
the boundary the temperature profile overshoots. 

choose to work in a 2D geometry to maximize the capability to get quantitative measure- 
ments at high Reynolds and Rayleigh numbers. Choosing a RT system presents the extra 
complexity of non-stationarity, allowing to probe also turbulent time scales. Our numerical 
simulations are performed using a recently proposed thermal Lattice Boltzmann Method [18] 
able to reproduce the Navier-Stokes equations for momentum, density and internal energy. 
Validation of the method can be found in [13]. Here we present numerical simulations up 
to 4096 x 10000 grid points. Table I provides details of our numerical experiment, obtained 
running on the QPACE supercomputer |19j . We limit our discussion here to the case of an 
ideal gas, with the equation of state P = pT, and at small Mach number. In this case, the 
main effect of stratification is limited to the presence of an adiabatic gradient affecting the 
evolution of temperature [H]. The equations are the following: (double indexes are summed): 




i = 



(1) 



1.06 




FIG. 2: RUN (B). Mean temperature profile T{z) after the RT evolution has stopped at the 
adiabatic slope (solid line). Bottom left panel: zoom of the overshooting region, with temperature 
inversion. Top right inset: heat flux profile u z 6{z) at three times before, during and after stopping: 
two regions develp with negative heat flux in correspondence of the temperature overshooting. 





L x L z 7 Ra m ax N con f L 7 t^r 


(A) 


4096 10000 -1-10" 5 8-10 9 18 10000 6.4 • 10 4 


(B) 


3072 7200 -4.2 • 10~ 5 3 • 10 9 11 2408 2.7 • 10 4 



TABLE I: Run (A): weak stratification. Run (B) strong stratification. Adiabatic gradient: 7 = 
-g/cp] c p = 2. Adiabatic length in grid units: L 7 = AT/I7I, AT = T down - T up , T m = (T down + 
T up )/2 and T up = 0.95, Td OW n = 1-05. Rayleigh number in presence of stratification is defined as [9]: 
Ra(t) = gL(t) 4 (AT/ L(t) + r y)/(uK)). Maximal value is obtained when L(t) = 3/4L 7 . Number of 
independent runs with random initial perturbation: N con j. Atwood number = At = AT/(2T m ) = 
0.05. Characteristic time scale, tRT = \/ L x / (gAt). 

where p is the deviation of pressure from the hydrostatic profile, p = P — Ph and d z Pu = 
—gpHi T m and p m are the mean temperature and density in the system, g is gravity and 
the adiabatic gradient is given by its ideal gas expression 7 = —g/c p with c p the specific 
heat. In this Boussinesq approximation j9] for stratified flows, momentum is forced only by 
temperature fluctuations 9 = T — T where we use the symbol (•) to indicate an average over 
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the horizontal statistically-homogeneous direction. The initial configuration is given by two 
8000 i , , , , , , 1 




4 6 8 t/t RT 12 14 



FIG. 3: Evolution of the mixing layer extension L{t) for both run (A) and (B). Notice the stop 
by stratification effects for the latter. Solid lines correspond to the closure prediction for the 
mixing length evolution obtained with the choice [1.45,4.5] (run A) and [0.3,6.5] (run B) for the 
parameters [be,bo Uz ] (see text). We also show 5 consecutive snapshots of the overshooting region 
highlighting the formation of a turbulent plume trying to enter the stable region and rejected back 
by gravitational forces. 



regions of cold (top) and hot (bottom) fluids prepared in the two half volumes of our cell 
(see right panel of Figjl]); turbulence is triggered by a small perturbation of the interface 
between them [20] (see [21] for a recent high resolution study of the classical RT system in 
3D). 

From Q, one easily derives the equations for the mean temperature profile, total kinetic 
energy, heat flux and temperature fluctuations: 

d t r + d z ^e = nd zz T 

\d t ¥ + \d z ffhT z + ^6(d z T - 7 ) = tM~e 

< _ (2) 

\d t u 2 + d z [u 2 u z + u z p] = -gdu z - t v 

d t W z + d z [6^ z + e^} = -g¥ + (3(z)^ z - e e , Uz 

V 
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where u 2 = u 2 x + it 2 ., (3(z) = (d z T — 7), and the dissipative terms are: ee >Uz = [y + n)di9diU z , 



eg = ndiOdiO, t v = vdiUjdiUj. The above set of equations is exact, except for boundary 



dissipative contributions as for example, nd z 6d z 6 which are irrelevant when k, v — V in 
absence of walls. From the second of ^ it is evident that temperature fluctuations are not 
forced anymore as soon as the mixing region develops a mean temperature profile of the 
order of the adiabatic gradient: 

d z T ~ 7. (3) 

As a consequence, from that time on, turbulence will decline. Given 7, we can identify the 
adiabatic length, L 7 = AT/I7I which corresponds to the prediction for the largest possible 
extension of the mixing layer during the RT evolution. In Fig. [3] we show the evolution of 
the mixing layer extension for two cases with weak (RUN A) and strong (RUN B) adiabatic 
gradient. Clearly, the extension of the mixing region stops to grow when L(t) ~ L 1 . In this 
paper we measure the mixing layer width L(t) as in [17J: 



L(t) =2 J dzB 



(4) 



with 0[f] = 2£; < £ < 1/2 and 9[f] = 2 (1 - £); 1/2 < £ < 1. The overshoot developing 
at the edge between turbulent and unmixed fluid is visible in Fig. |2j where we show both 
the nonlinear temperature inversion (inset bottom left) and the corresponding inversion in 
the heat flux (inset top right). This overshooting region is clearly a problem for any attempt 
to close the mean profile evolution with any sort of local eddy diffusivity: 

^J) = -K(z,t)d z T. 

Different models have been proposed in the literature for K(z, t), going from simple homoge- 
neous closure (K(z,t) oc L(t)\L(t)\) to Prandtl-like mixing-length closure [10l[TT] (K(z,t) oc 
L(t) 5 l 2 d z T) or following Spiegel's nonlinear approach [0] (K(z,t) oc L(t) 2 |9 z T| 1 / 2 ). All these 
attempts work well in the absence of overshooting and all of them suffer whenever an inver- 
sion in temperature and heat flux is observed (as in Fig. [2]), implying a negative effective 
eddy diffusivity. In order to overcome this difficulty, we keep exact the equations for the 
mean profile and close only the fluctuations at the second order moments in ([2]). Considering 
u z ~ u, we are left with six unknown to be defined: the three dissipative contributions on 
the rhs, and the three non-linear third order fluxes on the lhs. We close them adopting the 
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FIG. 4: Check of the closure for the case with weak stratification (run A). Circles: numerical 
data, solid lines: clousre. We plot the local heat flux vs the local temperature mean gradient, 9u z 
vs (d z T) x 10 5 , i.e. the effective diffusivity K(z,t) (). Error bars are evaluated averaging over 
the configurations specified in table (1). Result from the closure is given by the solid line with 
[b$,bg Uz ] = [1.45,4.5] and b Uz = 0.03. Inset: superposition of temperature profiles at three different 
times, t/tRT = 2,4,5.5 with the corresponding curves from the closure. 



simplest dimensionally-consistent local closure, for both fluxes and dissipative terms: 

9hT z = a e L\L\d zz ¥- e e = b e ¥/r(z, t) 

< (u 2 +p)u z = a Uz L\L\d zz u 2 ; e u = b Uz u 2 /r(z,t) 
0(u 2 z +p) = a dUz L\L\d zz 9u z ; e 6tUz = b 0Uz 6u z /T(z,t) 

where the typical time defining the dissipation rates is fixed by r(z,t) = \Jv? z /L{t). Some 
comments are in order. First, we notice that the closure is now local but on the second 
order moments, i.e. non-local for the evolution of mean profiles. Furthermore, out of the 
6 free parameters, 4 can be handled quite robustly, all the three coefficients in front of the 
closure for third order quantities are set 0(1) using a first order guess from the numerics, 
ag = 0.2; a Uz = 0.3; ag Uz = 0.8. Moreover, the free parameter defining the intensity of the 
kinetic dissipation, e V) is irrelevant in 2D (absence of direct energy cascade). The only 
two delicate free parameters are those defining the intensities of temperature and heat- 
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FIG. 5: Check of the closure for the case with strong stratification (run B). Circles: numerical 
data, solid line: closure. We plot the local heat flux vs the local temperature mean gradient, 0u z vs 
(d z T) x 10 5 , i.e. the effective diffusivity K(z,t) for two different times, t = 3tRT (bottom panel); 
t = QtfiT (top panel). Closure predictions ([be, be Uz ] = [0.35; 6.5]) are given by the solid lines. 
Inset bottom panel: matching between the temperature profile and the closure. Inset top panel: 
overshooting region around the top boundary layer. The two lines correspond to the closure using 
two different choices, [be,bg Uz ] = [0.25; 4.5]; [0.35,6.5]. 

flux dissipative terms, [be,be Uz ]- In order to get a good agreement with the numerics we 
need some fine tuning for them. It is important to notice that both dissipative terms will 
develop a non-trivial vertical profile, i.e. they are not given by a simple bulk homogeneous 
parameterization. 

In Fig. [3] we show that the closure is able to reproduce quantitatively the evolution of 
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L(t) for both cases of weakly stratified turbulence (A) and strongly stratified case (B). In 
Figs. (|4]|5]) we show the capability of the model to reproduce the heat flux vs. temperature 
profile behavior, providing a sort of aposteriori effective eddy diffusivity. As one can see, 
the closure is able to capture both the weak stratification case (Fig. [4]) and the strong 
stratification case (Fig. [5]). For the weak stratification case the aposteriori eddy diffusivity 
is very close to the one predicted by using the Prandtl mixing length theory as proposed in 
[TO] or to the one proposed by Spiegel's [9] (not shown). For the strong stratification case, 
our model is able to capture also the long time behavior, even after the evolution has come 
to a halt due to the adiabatic gradient, where the heat flux has completely inverted sign, 
see top panel of Fig. |5j In the inset of the same figure we also show the capability of the 
closure to reproduce the overshooting profile. As one can see the agreement is very good 
and the results are not very sensitive to the choice of the free parameters. 
We also remark that from eqs. (|2| one may derive a dimensional closure neglecting all terms 
except the one at large scales, and assuming that temperature fluctuations are dominated 
by the global temperature jump: 9 2 ~ AT. From the third equation we derive: 9u z ~ 
Lu 2 d z T and from the second: v? z ~ gL 2 d z T that -combined together- leads to Spiegel's 
closure: 9u z ~ L 2 (d z f)^ 2 d z f. 

In conclusions, we have performed state-of-the-art 2D numerical simulations using a novel 
LBM for turbulent convection driven by a Rayleigh- Taylor instability in weakly and strongly 
stratified atmospheres. For the strongly stratified case, we are able to resolve the overshoot- 
ing region with up to 800 grid points, something impossible to achieve with 3D direct nu- 
merical simulations because of lack of computing power. We have presented a second-order 
closure to describe the evolution of mean fields able to capture both bulk properties and 
the overshooting observed at the boundary between the stable and unstable regions. The 
closure is local in terms of fluctuations while keeping the evolution of the mean temperature 
profile exact. In order to apply the same closure to 3D cases one needs to take into account 
some possible non-trivial effects induced by the kinetic energy dissipation modeling, due to 
the presence of a forward energy cascade. As a result, also the parameter b Uz will become a 
relevant input in the model. 

We acknowledge useful discussions with G. Boffetta and A. Wirth. We acknowledge 
access to QPACE and eQPACE during the bring-up phase of these systems. Parts of the 
simulations were also performed on CASPUR under HPC Grant 2009 and 2010. 
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